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Abstract 

Anthropogenic activities frequently result in both rapidly changing environments 
and translocation of species from their native ranges (i.e., biological invasions). 
Empirical studies suggest that many factors associated with these changes can lead 
to complex genetic patterns, particularly among invasive populations. However, 
genetic complexities and factors responsible for them remain uncharacterized in 
many cases. Here, we explore these issues in the vase tunicate Ciona intestinalis 
(Ascidiacea: Enterogona: Cionidae), a model species complex, of which spA and spB 
are rapidly spreading worldwide. We intensively sampled 26 sites (N = 873) from 
both coasts of North America, and performed phylogenetic and population genet- 
ics analyses based on one mitochondrial fragment (cytochrome c oxidase subunit 3 
-NADH dehydrogenase subunit I, COX3-ND1) and eight nuclear microsatellites. 
Our analyses revealed extremely complex genetic patterns in both species on both 
coasts. We detected a contrasting pattern based on the mitochondrial marker: two 
major genetic groups in C. intestinalis spA on the west coast versus no significant 
geographic structure in C. intestinalis spB on the east coast. For both species, geo- 
graphically distant populations often showed high microsatellite-based genetic 
affinities whereas neighboring ones often did not. In addition, mitochondrial and 
nuclear markers provided largely inconsistent genetic patterns. Multiple factors, 
including random genetic drift associated with demographic changes, rapid selec- 
tion due to strong local adaptation, and varying propensity for human-mediated 
propagule dispersal could be responsible for the observed genetic complexities. 
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Introduction 

Identifying factors responsible for microevolutionary pro- 
cesses over contemporary timescales is important for under- 



standing how species evolve in response to rapidly changing 
environments associated with recent anthropogenic activi- 
ties (Carroll 2008; Roux and Wieczorek 2008). Biological 
invasions provide an opportunity to answer fundamental 
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questions regarding microevolution, especially those involv- 
ing rapid range expansion or habitat transitions (Lee 2002; 
Suarez and Tsutsui 2008; Winkler et al. 2008; Bock et al. 
2012). Processes associated with biological invasions, such 
as transfers of individuals to dramatically different environ- 
ments and sudden changes in population demography, can 
lead to exceedingly rapid evolutionary changes (Reznick and 
Ghalambor 2001; Lee 2002; Prentis et al. 2008). However, 
understanding evolutionary processes that may accompany 
biological invasions can be challenging because of genetic 
complexities generated by biological invasions (Lee 2002). 
Despite recent progress on evolutionary aspects of biological 
invasions, the role and importance of the microevolutionary 
factors remain poorly understood (see reviews by Lee 2002; 
Prentis et al. 2008). 

To better understand these evolutionary mechanisms, it 
is necessary to identify potential generalities and inconsis- 
tencies within and among multiple invasive species across 
different geographic scales, and then identify the evolution- 
ary/ecological factors associated with common and unique 
genetic patterns. For aquatic species, many factors including 
life-history characteristics (e.g., sexual vs. asexual reproduc- 
tion, planktotrophic vs. direct development, etc.), environ- 
mental heterogeneity (e.g., different temperature, salinity, 
etc.), anthropogenic activities (e.g., habitat fragmentation, 
pollution, etc.), coupled with varying propensity for human- 
mediated propagule dispersal, can lead to extremely complex 
genetic patterns (Lodge 1993; Kolar and Lodge 2001; Sakai 
et al. 2001; Lee and Boulding 2009; Dupont et al. 2010; Bock 
et al. 2011). Recent empirical studies have revealed complex 
genetic patterns at different geographic scales in aquatic in- 
vasive species (see e.g., Ben-Shlomo et al. 2006; Dupont et al. 
2009; Tepolt et al. 2009; Bock et al. 201 1). Given that numer- 
ous nonindependent factors may be involved in determining 
genetic patterns of aquatic invasive species, the choice of an 
appropriate model system seems critically important. 

The vase tunicate, Ciona intestinalis, is a well-known model 
species for developmental and evolutionary biology (see re- 
view by Procaccini et al. 20 1 1 ) . The recently published whole 
genome (Dehal et al. 2002), coupled with detailed phylo- 
genetic and population genetics analyses (Caputi et al. 2007; 
Sordino et al. 2008; Nydam and Harrison 2010; Zhan et al. 
2010), has expanded C. intestinalis as a model for evolution- 
ary studies associated with biological invasions. Phylogenetic 
analyses revealed that C. intestinalis is a species complex con- 
sisting of at least four morphologically cryptic but genetically 
distinct species (spA, spB, spC, spD; Zhan et al. 2010). Inter- 
estingly, the two highly invasive species, spA and spB, inhabit 
largely disjoint geographic regions worldwide, while the other 
two, spC and spD, remain restricted to their putative native 
ranges in the Mediterranean Sea and Black Sea, respectively 
(Caputi et al. 2007; Nydam and Harrison 2007; Zhan et al. 
2010). Population genetics analyses further illustrated rel- 
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atively low population differentiation and high population 
connectivity for the two highly invasive species at both re- 
gional and continental scales (Zhan et al. 2010). Abundant 
genetic information generated by intensive developmental 
and evolutionary studies, as well as increasing knowledge re- 
garding phylogenetics and population biology, indicates that 
the C. intestinalis species complex is a promising system for 
comparing genetic patterns in closely related invasive species 
during range expansions. 

During the past century, the two highly invasive species 
of C. intestinalis complex, spA and spB, have successfully 
invaded coastal marine habitats throughout the temperate 
zone (Kott 1952; Monniot and Monniot 1994; Lambert and 
Lambert 2003; McDonald 2004). In invaded areas, these two 
species rapidly cover available substratum (Fig. 1) and ex- 
clude native species, often causing considerable economic 
damage and ecological changes (Lambert and Lambert 2003; 
Ramsay et al. 2008). In North America, C. intestinalis col- 
onized the west coast either in the early 1910s (Huntsman 
1912) or in the 1930s (Lambert and Lambert 1998). Un- 
certainty in its invasion history is largely due to the lack 
of distinct morphological attributes between C. intestinalis 
and its congener, C. savignyi (Lambert 2003). Owing to the 
absence of records on the east coast, the invasion history 
remains uncertain along almost all of this coast, with the ex- 
ception of bays in Prince Edward Island, which were newly 
invaded in 2004 (Ramsay et al. 2008). So far, C. intestinalis 
has successfully colonized the Pacific coast of North America 
from British Columbia, Canada, to southern California, 
USA, and the Atlantic coast from Newfoundland, Canada, 
to Florida, USA (Therriault and Herborg 2008; de Oliveira 
Marins etal. 2009). 

Here, we use phylogenetic and population genetics ap- 
proaches based on mitochondrial DNA (mtDNA) sequences 
and nuclear microsatellites to characterize genetic patterns of 
26 populations of C. intestinalis sampled from both coasts of 
North America. We aim to clarify population genetic struc- 
ture of the two highly invasive species, C. intestinalis spA and 
spB on both coasts of North America. Based on these analy- 
ses, we discuss the potential factors driving complex genetic 
patterns, especially among closely related invasive species. 

Materials and Methods 

Sampling and species identification 

Ciona intestinalis specimens were collected from the invaded 
ranges on both coasts of North America. On the west coast, 
sampling was conducted along the coast of California be- 
tween July and November 2006, spanning from Tomales 
Bay to San Diego estuary, while on the east coast sampling 
was performed between October 2007 and March 2009 in 
Prince Edward Island, Nova Scotia, and Connecticut (Fig. 2). 
In total, we sampled 26 populations (N = 873; Fig. 2; 
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Figure 1. The vase tunicate, Ciona intestinalis, fouled with high density on hard substratum in the Port of Point Tupper (left upper), on sampling 
plates (before and after sampling, left lower), and on ropes for mussel aquaculture in Prince Edward Island (right). Photos by Anai's Lacoursiere-Roussel 
and Samuel Collin. 




Figure 2. Sampling sites and distribution of mitochondrial cytochrome c oxidase subunit 3-NADH dehydrogenase subunit 1 (COX3-ND1) haplotypes 
for Ciona intestinalis spA on the west coast and spB on the east coast of North America. Site IDs as per Table 1 . Pie charts indicate the proportion of 
haplotype groups observed at each site, with different colors corresponding to different haplotypes. 



Table 1), including 11 from the west coast (N = 380) and 
15 (N = 493; including populations BR, MR, PT, CT, LT, and 
GT from Zhan et al. 2010) from the east coast. All specimens 
were preserved in 95% ethanol prior to genetic analyses. 



Owing to the lack of reliable diagnostic morpho- 
logical characters among species in C. intestinalis complex, 
all collected specimens were identified to species level us- 
ing one mtDNA fragment, cytochrome c oxidase subunit 
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3-NADH dehydrogenase subunit 1 (COX3-ND1; Zhan et al. 
2010). 

mtDNA amplification and sequencing 

Total genomic DNA was extracted from approximately 50 mg 
of tissue according to the proteinase K method (Waters et al. 
2000). The mitochondrial COX3-ND1 fragment was ampli- 
fied using the primers TX3F and TN1R (Iannelli et al. 2007). 
PCR amplification was performed according to the proto- 
col described by Iannelli et al. (2007). Sequencing reactions 
were performed using the forward primer (TX3F), BigDye 
Terminator v3.1 sequencing chemistry, and an ABI 3130XL 
automated sequencer (PE Applied Biosystems, Foster City, 
CA). Sequences that contained ambiguous nucleotides were 
subsequently sequenced using the reverse primer (TN1R). 

mtDNA analysis 

Mitochondrial sequences were aligned and edited us- 
ing CodonCode Aligner v2.0.6 (CodonCode Corporation, 
Dedham, MA). Bayesian inference (BI) and neighbor- 
joining (NJ) phylogenetic analyses were conducted using 
C. savignyi as outgroup (Nydam and Harrison 2007; Zhan 
et al. 2010). The best DNA substitution model for mtDNA 
haplotypes was determined using MODELTEST v3.7 (Posada 
and Crandall 1998) with Akaike Information Criterion. The 
BI analysis was conducted using MRBAYES v3.2 (Ronquist 
and Huelsenbeck 2003). Trees were sampled every 100 gen- 
erations for two million generations, and the first 25% of 
all the trees sampled were discarded as burn-in. NJ phylo- 
genetic analyses were performed using MEGA v4.0 (Tamura 
et al. 2007) based on nucleotide distances corrected using the 
Tamura-Nei model (Tamura and Nei 1993). Clade support 
was estimated using bootstrap analysis with 1000 replicates. 
Relationships between mtDNA haplotypes were further ex- 
amined using a statistical parsimony haplotype network gen- 
erated at the 95% connection limit with TCS vl.2 1 (Clement 
et al. 2000). 

The number of haplotypes («), haplotype diversity (h), 
and nucleotide diversity (ji ) were assessed using DNASP v5 
(Rozas et al. 2003). Population pairwise <t> S T was calculated 
using the Tamura and Nei (TrN) substitution model and 
10,000 permutations in ARLEQUIN v3.1 (Excoffier et al. 
2005), with levels of significance adjusted using sequential 
Bonferroni corrections (Rice 1989). To assess hierarchical 
population genetic structure, we conducted a hierarchical 
analysis of molecular variance (AM OVA; Excoffier et al. 1 992 ) 
based on 10,000 random permutations using ARLEQUIN. 
AMOVA was performed separately for both coasts (i.e., two 
species, spA on the west coast and spB on the east coast, 
see Results section). Populations were grouped according to 
different geographical regions on both coasts, two groups 
on the west coast: northern California (TB, SF, MO) and 



southern California (SB, CI, PH, LA, NB, OE, MI, SD; 
Table 1, Fig. 2); three groups on the east coast: Prince Edward 
Island (CR, BR, MR), Nova Scotia (SN, PT, HF, CT, MA, 
MB, ST, LU, SH, LT, YM), and Connecticut (GT; Table 1; 
Fig. 2). Molecular variance was partitioned into three lev- 
els: between regions, among populations within regions, and 
within populations. To test the correlation between genetic 
distance [ < 1 > st/(1- < 1 > st)] and geographic distance, we per- 
formed isolation by distance (IBD) analysis for spA and spB 
using a Mantel test with 10,000 permutations implemented in 
GENEPOP v3.4 (Raymond and Rousset 1995). Geographical 
distances were calculated as the minimum coastline distances 
between sampling locations using Google Earth. 

Microsatellite DNA genotyping and analysis 

All individuals were genotyped for eight polymorphic 
microsatellite markers (Zhan et al. 2010). PCR amplifica- 
tion and genotyping were performed according to the pro- 
tocol described by Zhan et al. (2010). Microsatellite genetic 
diversity was evaluated by the number of alleles (A), allelic 
frequency (F), allelic richness (A r ), as well as the observed 
heterozygosity (H 0 ) and expected heterozygosity (H E ) using 
Fstat v2.9.3.2 (Goudet 2001). Additionally, a nonparametric 
test (Mann-Whitney U test), implemented in STATISTICA 
v6 (StatSoft, Inc., Tulsa, OK), was used to test the difference in 
allelic richness (A r ) and expected heterozygosity (H E ) among 
populations. Markov chain method (Guo and Thompson 
1992) was employed to estimate the probability of signifi- 
cant deviation from Hardy-Weinberg equilibrium (HWE) 
using GENEPOP. Significance criteria were adjusted for the 
number of simultaneous tests using sequential Bonferroni 
corrections. 

Population differentiation was determined by F S j (Weir 
and Cockerham 1984) for all population pairs using F S tat- 
A total of 10,000 permutations were performed, and signif- 
icance levels were adjusted using sequential Bonferroni cor- 
rections. In addition, population structure was determined 
by conducting a three-dimensional factorial correspondence 
analysis (3D-FCA) in GENETIX v4.05 (Belkhir et al. 2004) 
and a Bayesian, Markov Chain Monte Carlo (MCMC) based 
approach implemented in STRUCTURE v2.1 (Pritchard et al. 
2000). For the STRUCTURE analysis, we assessed likelihoods 
for models with the number of clusters ranging from K = 1 to 
the total number of sampling sites (11 forspAand 15for spB). 
Ten independent runs were performed for each specified re- 
value, and for each run, 1,000,000 generations were used after 
discarding 100,000 generations as burn-in. The measure of 
AK (Evanno et al. 2005) was used to infer the number of 
biologically relevant clusters. To test for a pattern of IBD, we 
used a Mantel procedure with 10,000 permutations imple- 
mented in GENEPOP to assess the dependence between the 
genetic distances [fsT/(l-fsT)] and geographic distances in 
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each species. A hierarchical AMOVA was performed in AR- 
LEQUIN based on microsatellite genetic distances between 
populations, and partitioning variance as for mtDNA. 

Results 
mtDNA analyses 

A total of 663 individuals derived from 26 populations from 
both coasts were successfully sequenced for the COX3-ND1 
fragment (Table 1). Species identification based on phylo- 
genetic analysis of this fragment recovered two well- 
supported phylogroups, corresponding to spA and spB 
(Fig. 3A). The mean interclade divergence was 15.9%, much 
higher than the intraclade divergences (spA: 0.5%; spB: 
0.8%). All individuals from west coast were identified as spA, 
while all individuals from east coast corresponded to spB 
(Table 1; Fig. 3A). 



In total, we identified 19 haplo types (GenBank accession 
numbers: HM036361, HM036366, JQ396394-JQ396310) 
for spA and 27 haplotypes (GenBank accession num- 
bers: HM036368, HM036370, HM036375-HM036381, 
HM036383-HM0385, JQ396311-JQ396325) for spB 
(Table 1; Fig. 2). The number of haplotypes per population 
ranged from two to eight for spA and three to eight for 
spB, while the haplotype diversity (h) varied from 0.429 to 
0.861 for spA, and 0.191 to 0.821 for spB (Table 1). The 
lowest genetic diversity (h = 0.429) in spA was found in the 
population MO (Monterey Bay), sampled from the open 
sea shore on the west coast, while the highest was detected 
in major ports including SF (San Francisco, h = 0.861), PH 
(Port Hueneme, h = 0.752), and SD (San Diego, h = 0.714; 
Table 1). Similarly, some populations of spB sampled from 
ports showed relatively high genetic diversity (e.g., h = 0.577 
for population Sydney; Table 1), but some did not (e.g., 
h = 0.267 for population Point Tupper; Table 1). Two of the 



Figure 3. Bayesian inference (Bl) tree (A) 
based on mitochondrial cytochrome c oxidase 
subunit 3-NADH dehydrogenase subunit 1 
(COX3-ND1) haplotypes, and haplotype 
networks generated with TCS for Ciona 
intestinalis spA (B) on the west coast and 
C. intestinalis spB (C) on the east coast of 
North America. Haplotype and population 
names as per Table 1 . For the Bl tree, posterior 
probabilities for Bayesian inferences (in 
percentage, above branch) and bootstrap 
values (below branch) for neighbor-joining 
reconstruction are shown at major nodes. For 
the TCS network, sampled haplotypes are 
indicated by circles and missing or unsampled 
haplotypes are indicated by solid black dots. 
Haplotypes are color coded according to 
different geographical regions. 
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most recently established populations, BR (Brudenell River; 
reported in 2005, C. D. Mills, personal communication) and 
MR (Murray River; detected in 2006, C. D. Mills, personal 
communication), sampled from aquaculture facilities in 
Prince Edward Island, exhibited the lowest haplotype 
diversity (h = 0.191; Table 1). By contrast, one recently 
established neighboring population CR (Cardigan River; 
reported in 2006, C. D. Mills, personal communication) 
showed a similar level of genetic diversity (h = 0.556) to that 
in other populations (Table 1). 

The shallow phylogeny within each species was confirmed 
by TCS haplotype networks. Haplotypes were connected to 
each other by only several mutation steps (Fig. 3B and C) . The 
network for spB showed a star-shaped pattern with one dom- 
inant haplotype at the center (Fig. 3C), whereas we detected 
two dominant haplotypes at the center for spA (Fig. 3B). De- 
spite the shallow phylogeny and close relationships among 
haplotypes, we detected geographic structure for spA along 
the west coast. Generally, populations located at two sides of 
Point Conception belonged to two geographic groups, that 
is, southern and northern California (Fig. 2). The haplo- 
type composition and haplotype frequencies were different 
in these two geographic groups (Table 1; Fig. 2). For example, 
the most dominant haplotype in southern California, Ca7, 
was not detected in any locations in northern California. Sim- 
ilarly, haplotypes with relatively high frequencies in north- 
ern California, such as Ca20, were not observed in southern 
California (Fig. 2). The difference between two groups was 
further confirmed by both pairwise <J> S t values (Table 2 A) and 
AMOVA (Table 3). Comparisons between populations sam- 
pled from northern and southern California revealed high 
and significant O st values in more than half of the popu- 
lation pairs (14 out of 27; Table 2A), while low and non- 
significant values were detected for population pairs within 
each region (Table 2A). AMOVA attributed a statistically sig- 
nificant proportion of genetic variance (14.53%, P < 0.005) 
to the among-group component (i.e., between two groups, 
Table 3). 

By contrast, we did not observe significant geographic 
structure at the mtDNA locus for spB along the east 
coast. Generally, all populations shared the same dom- 
inant haplotype, Cb4 (Table 1; Fig. 2). Statistical anal- 
ysis of pairwise <J> S t values revealed that only eight of 
120 pairs were significantly different after sequential Bon- 
ferroni correction (Table 2B). AMOVA attributed only a 
small proportion of genetic variance (3.42%) to the among- 
group component, although it was statistically significant 
(Table 3). 

Mantel tests based on the COX3-ND1 fragment revealed 
no correlation between genetic and geographical distances for 
either species on both coasts (Fig. 4), suggesting no pattern 
oflBD. 
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Microsatellite DNA analyses 

A total of 709 individuals from 26 populations, including 301 
individuals from 1 1 populations for spA and 408 individu- 
als from 15 populations for spB, were successfully genotyped 
at eight microsatellite loci. In total, we identified 183 and 
239 alleles for spA and spB, respectively. Allelic richness (A r ) 
ranged from 3.4 to 7.0 (mean 5.3) for spA and from 3.7 to 9.3 
(mean 6.8) forspB (Table 1). Generally, high expected hetero- 
zygosity (H E ) values were observed for all populations, rang- 
ing from 0.563 to 0.703 for spA and from 0.510 to 0.875 
for spB (Table 1; Table SI). For spA, the population MO 
(Monterey Bay) had the lowest genetic diversity, all consis- 
tent with mtDNA results (Table 1). Allelic richness in this 
population was significantly lower than three populations 
(LA, OE, SD) (Mann-Whitney U test, P < 0.013; Table S2). 
We did not detect statistical difference in expected hetero- 
zygosity between any population pairs (P > 0.05; Table S2). 
Similarly to mtDNA, populations sampled from major ports 
showed relatively high genetic diversity (Table 1). 

Interestingly, nuclear and mtDNA markers often provided 
inconsistent results. For example, for spA population PH 
(Port Hueneme) had the highest number of mtDNA haplo- 
types (n = 8; Table 1) but a relatively low number of alleles 
at microsatellite loci (mean A = 4.1; Table 1). For spB, the 
lowest microsatellite genetic variation was identified in four 
populations including MA (Martin's River), MB (Mahone 
Bay), ST (Stone Hurst), and LU (Lunenburg), which were 
sampled from southeastern Nova Scotia (Tables 1 and SI). 
Both allelic richness and expected heterozygosity were signifi- 
cantly lower than those in other populations (Mann- Whitney 
U test; Table S2 ) . However, two of these four populations, MA 
and ST, exhibited high mitochondrial diversity (Table 1). 
Compared to the other populations, the populations sam- 
pled from major ports such as HF (Halifax) did not show 
higher genetic diversity, and recently established populations 
in Prince Edward Island did not exhibit significantly lower 
genetic diversity (Table 1). Deviations from HWE were ob- 
served at multiple loci and sampling locations for both species 
(Table SI). All of the deviant cases showed significant het- 
erozygote deficiency (P < 0.001; Table SI). As suggested by 
Zhan et al. (2010), recurrent inbreeding and Wahlund effects 
may be responsible for these deviations. 

Overall, we detected a relatively high level of genetic dif- 
ferentiation in both species when analyzing microsatellite 
data. For spA, pairwise F S t values ranged from zero to 0.123 
with an average of 0.059. In total, 36 of 55 (65.5%) popu- 
lation pairs were significantly different after Bonferroni cor- 
rection (Table 2A). Interestingly, the highest genetic differen- 
tiation (_Fst = 0.123) was detected between two neighboring 
populations, CI (Channel Islands) and PH (Port Hueneme; 
Fig. 2; Table 2A). For spB, pairwise F S j values varied from 
0.005 to 0.191, with an average of 0.088. In total, 75 of 105 
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Table 2. Estimates of population genetic differentiation in two highly invasive vase tunicates: Ciona intestinalis spA on the west coast (A) and spB on 
the east coast (B) of North America. Above diagonal: pairwise <i> ST based on cytochrome c oxidase subunit 3-NADH dehydrogenase subunit 1 region 
(COX3-ND1); Below diagonal: pairwise F ST based on microsatellite markers. Bold numbers indicate statistical significance after sequential Bonferroni 
corrections. Negative values were converted into zero. 
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0.265 
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(71.4%) comparisons of population pairs remained signifi- 
cantly different after Bonferroni correction (Table 2B). The 
highest P S t values were detected for comparisons between 
four populations (MA, MB, ST, LU) and all other sites sur- 
veyed, with an average of pairwise P ST values greater than 0. 1 
(Table 2B). 

The two different approaches used to infer population 
structure, 3D-FCA and Bayesian clustering, provided largely 
consistent results in both species. For spA on the west 
coast, both analyses grouped populations into two clusters 
based on geographic regions (i.e., northern and southern 
California), which is inconsistent with mtDNA results ob- 
tained from pairwise <J> S t (Table 2A) and AMOVA (Table 3). 
Bayesian clustering analysis suggests a two-cluster model 
(K = 2) as the most parsimonious possibility (Fig. 5A). 
Most individuals from five populations (TB, SF, CI, LA, MI) 
were assigned into one cluster regardless of geographic ori- 
gin, whereas most individuals from three populations (SD, 
SB, OE) collected from southern California were assigned to 
the other cluster (Fig. 5A). However, individuals from the 



remaining two populations, PH (Port Hueneme) and MO 
(Monterey Bay), do not appear to assign consistently to a 
single cluster (Fig. 5A). This pattern was confirmed by 3D- 
FCA. Component 1, explaining 22.95% of genetic variance, 
separated the five (TB, SF, CI, LA, MI) from the three popula- 
tions (SD, SB, OE), with the remaining two (PH, MO) in the 
middle (Fig. 6A). For spB on the east coast, the two analyses 
consistently supported a two-cluster model and clearly sepa- 
rated four populations (MA, MB, ST, LU; Figs. 5B and 6B), all 
consistent with microsatellite diversity results (Table 1) and 
pairwise Fst values (Table 2B). Since low genetic diversity 
in these four populations can inflate genetic differentiation, 
3D-FCA and Bayesian clustering were rerun after removing 
these four populations. Bayesian clustering analysis still sup- 
ported a two-cluster model, showing a high level of admixture 
and mixture of individuals from different geographical ori- 
gins (Fig. 5C). This finding is consistent with the 3D-FCA 
results, showing no obvious population cluster (Fig. 6C). 
Similar to spA, the two clusters showed discontinuous distri- 
butions along the east coast. Neighboring populations were 
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Table 3. Results of the analysis of molecular variance (AMOVA) for two highly invasive species, Ciona intestinalis spA on the west coast and spB on 
the east coast of North America. Populations were grouped according to different geographic regions on both coasts, two groups on the west coast: 
northern California (TB, SF, MO) and southern California (SB, CI, PH, LA, NB, OE, Ml, SD); three groups on the east coast: Prince Edward Island (CR, 
BR, MR), Nova Scotia (SN, PT, HF, CT, MA, MB, ST, LU, SH, LT, YM), and Connecticut (GT). 
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Figure 4. Correlation between geographical distance in kilometers (x-axis) and genetic distance given as F ST /(1 - F ST ) for microsatellite data and 
ct> 5T /( 1 - 4> ST ) for mitochondrial data (y-axis) in Ciona intestinalis spA on the west coast (A) and spB on the east coast (B) of North America. Open and 
solid cycles represent data derived from mitochondrial and microsatellite markers, respectively. 
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Figure 5. Bayesian inference (Bl) of population structure for Ciona intestinalis spA (A), spB (B), and data subset of spB (C) based on microsatellite 
markers. For Bayesian clustering analysis, each genotype is represented by a thin vertical line, with proportional membership in different clusters 
indicated by colors. Bold vertical lines separate collection sites, with site IDs indicated below the plot. 
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Figure 6. Three-dimensional factorial correspondence analysis (3D-FCA) of Ciona intestinalis spA (A), spB (B), and data subset of spB (C) after removing 
four populations (MA, MB, ST, and LU) based on eight polymorphic microsatellites. 



not assigned to a single cluster (e.g., CT and SH), whereas 
geographically isolated populations were frequently assigned 
to the same clusters (e.g., GT and CR). In addition, Mantel 
tests based on microsatellite data showed no indication of 
IBD in both species (Fig. 4). 

Discussion 

Our genetic surveys based on both population genetics and 
phylogenetic analyses show complex genetic patterns in two 
closely related invasive species, C. intestinalis spA and spB. 
A contrasting geographic pattern was detected in the two 
species based on the mitochondrial marker: two major ge- 
netic groups in spA on the west coast versus no signifi- 
cant geographic structure in spB on the east coast (Fig. 2; 
Table 2). For both species on both coasts, multiple anal- 
yses showed that geographically distant populations were 
often genetically similar, whereas neighboring populations 
were sometimes not (Figs. 5 and 6). In addition, we ob- 
served inconsistent genetic patterns when using mitochon- 
drial versus nuclear markers. All these genetic complexities 



suggest that multiple factors, coupled with their interactions, 
could be responsible for the genetic patterns observed in 
Ciona. 

Geographic structure based on mtDNA in 
spA on the west coast 

Genetic analyses based on mtDNA revealed contrasting pat- 
terns between spA on the west coast and spB on the east 
coast: two major geographic groups in spA versus no signif- 
icant geographic structure in spB (Fig. 2). For spA on the 
west coast, the genetic patterns can result from three possible 
processes: ( 1 ) separate introductions from genetically distinct 
sources; (2) random genetic drift associated with population 
demography; and/or (3) rapid selection due to strong local 
adaptation. 

Given extremely heavy shipping to the west coast of 
North America (see e.g., Sylvester et al. 2011), it is possi- 
ble that separate introductions may have occurred in both 
Ciona and other ascidians (see e.g., Stoner et al. 2002; Zhan 
et al. 2010; Bocket al. 2011). However, analyses based on mi- 
crosatellites suggest that this process is unlikely responsible 
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for defining the two geographic groups. Two tests, including 
3D-FCA and Bayesian clustering, did not separate north- 
ern from southern populations (Figs. 5A and 6A). In addi- 
tion, microsatellite analyses such as pairwise F S j values sug- 
gest a high level of gene flow between several northern and 
southern populations. For example, there is no detectable 
microsatellite-based genetic differentiation between popula- 
tions TB (Tomales Bay) and NB (Newport Bay), located in 
two different geographic groups as defined by mtDNA (Fig. 2; 
Table 2). Additional evidence for supporting this argument 
comes from little genetic structure based on global sampling 
and population genetics and phylogeographic analyses us- 
ing both nuclear and mtDNA markers (Zhan et al. 2010). 
Altogether, available evidence suggests that separate intro- 
ductions are unlikely responsible for such a genetic pattern 
in spA on the west coast. 

Genetic drift can create complex genetic patterns by ran- 
domly changing allele frequencies and/or eliminating low- 
frequency alleles, especially when effective population sizes 
(N e ) decrease. Moreover, compared to nuclear markers, 
mtDNA has smaller effective population size (N emtDNA ~ 1/4 
JVenuciear)' likely leading to stronger signatures of genetic drift 
on mtDNA. For broadcast spawning marine invasive species 
including Ciona, two possible processes — sweepstakes re- 
productive success and dramatic demographic changes dur- 
ing translocation — can lead to strong genetic differentiation 
between populations. Marine species with broadcast spawn- 
ing exhibit sweepstakes reproductive success, in which a sub- 
set of the population contributes the majority of offspring 
to subsequent generations (Hedgecock et al. 2007). This can 
sharply decrease effective population size, resulting in spatial 
and/or temporary genetic structuring (see e.g., Li and Hedge- 
cock 1998, Hedgecock et al. 2007; Christie et al. 2010). The 
existence of population-specific haplotypes/alleles (Fig. 2; 
Table SI) suggests that sweepstakes reproductive success 
might occur in Ciona. In addition, population bottlenecks 
may have occurred during the translocation and settlement 
on the west coast. Random genetic drift due to severe bottle- 
necks, especially for mtDNA that is more sensitive to effective 
population size changes, may have created divergent genetic 
groups at initial stages of invasions when limited numbers 
of propagules were transplanted, and then these divergent 
genetic groups were transported to other locations via sec- 
ondary spread. 

However, genetic drift alone seems inadequate to explain 
the genetic pattern on the west coast, given the incon- 
sistency between unique dominant haplotypes to specific 
geographic regions (Table 1; Fig. 2) and high gene flow 
between these regions suggested by microsatellite analyses 
(Table 2A; Figs. 5A and 6A). This mitochondrial-nuclear 
discrepancy (i.e., high mtDNA differentiation vs. no or 
very low genetic differentiation at nuclear loci) has been 



detected in several other marine species including those 
with relatively high self-dispersal capacity (e.g., a marine 
goby Pomatoschistus minutus; Larmuseau et al. 2010). Se- 
lection and sex-biased dispersal are generally adopted as 
major explanations for this discrepancy (see Larmuseau 
et al. 2010 and references therein). However, given that 
Ciona is hermaphroditic, sex-biased dispersal cannot be 
a valid explanation. Considering a high level of ge- 
netic exchange among populations suggested by micro- 
satellite analyses (Table 2A; Figs. 5A and 6A) and possible 
population connection mediated by frequent shipping, such 
differences at mtDNA between northern and southern Cali- 
fornia may be maintained by rapid selection associated with 
strong local adaptation. The two geographic groups detected 
on the west coast are located on two sides of a well-known 
oceanographic and biogeographic boundary, Point Concep- 
tion (PC), where multiple species exhibit phylogeographic 
breaks (see e.g., Burton 1998; Bernardi 2000; Dawson et al. 
2006; Earl et al. 2010). PC is marked by strong discontinu- 
ities in water temperature, salinity, dissolved oxygen, as well 
as hydrography (Burton 1998; Dawson et al. 2001). Strong 
upwelling and high wave exposure dominate the north of PC, 
while weak and seasonal upwelling and warmer water tem- 
peratures are characteristics in the south (Blanchette et al. 
2007). These oceanographic conditions have important ef- 
fects on benthic communities, influencing both growth rate 
and recruitment dynamics (see e.g., Connolly et al. 2001; 
Blanchette et al. 2007). Therefore, Ciona populations might 
locally adapt to these two dramatically distinct environments. 
Compared to neutral nuclear makers such as microsatellites, 
selection events are much easier to capture through mtDNA, 
mainly because the entire mitochondrial genome is one link- 
age group without recombination. A selective sweep due to 
any selection pressure on any mtDNA gene can easily fix the 
haplotypes having higher fitness. Indirect and direct selec- 
tions on mtDNA have been demonstrated in many species 
(e.g., reviews by Ballard and Rand 2005; Hurst and Jiggins 
2005), although genetic mechanisms for such selections on 
mtDNA remain poorly understood. Consequently, the ob- 
served pattern on the west coast in spA likely resulted from 
the interaction between random genetic drift associated with 
demographic change and rapid selection due to strong local 
adaptation. 

Additional evidence suggests that the interaction between 
genetic drift and selection could be responsible for such a geo- 
graphic pattern. When our data are combined with the global 
data from Zhan et al. (2010), we observed different dominant, 
and sometimes unique, haplotypes on different continents, 
although only several mutation steps were observed between 
these haplotypes (Fig. SI). This pattern was found not only 
in spA but also in spB. For instance, all haplotypes of spB in 
Asia are unique (Fig. SI). 
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Complex intraspecific genetic patterns in 
both species 

Multiple analyses show that intraspecific genetic patterns var- 
ied significantly across the sampling regions analyzed. High 
genetic differentiation (Fst values) was observed not only be- 
tween geographically distant populations, but also between 
adjacent ones (Table 2). IBD analyses did not reveal any sig- 
nificant correlation between genetic differentiation and geo- 
graphical proximity for either species based on both types of 
markers (Fig. 4). Both Bayesian clustering and 3D-FCA occa- 
sionally assigned geographically distant populations into the 
same clusters (Figs. 5 and 6). 

Recent empirical studies have revealed that human- 
mediated dispersal may leave strong signatures on population 
genetic structures in numerous invasive species, including as- 
cidians (see e.g., Dupont et al. 2010; Zhan et al. 2010; Bock 
et al. 2011). Human-mediated gene flow can both generate 
complex genetic structuring at regional and local scales (see 
e.g., Darling and Folino-Rorem 2009; Dupont et al. 2010; 
Bock et al. 2011), and affect population genetic patterns 
among continents (see, e.g., Zhan et al. 2010). Generally, 
aquaculture transfer and shipping are considered as the lead- 
ing pathways for human-mediated dispersal of ascidians in 
North America (Dijkstra et al. 2007; Zhan et al. 2010; Bock 
et al. 2011). Our genetic analyses detected apparent signa- 
tures of these two pathways on population genetic structure 
in Ciona. The major pathway responsible for introductions of 
Ciona to Prince Edward Island is aquaculture trade between 
Prince Edward Island and Nova Scotia (Dijkstra et al. 2007). 
Our genetic analysis is concordant with this view. For ex- 
ample, the population CT (Chester), sampled from Chester 
basin supporting aquaculture in Nova Scotia, showed no sig- 
nificant genetic differentiation from the three populations 
(CR, BR, MR) sampled from aquaculture facilities in Prince 
Edward Island (Table 2). Genetic analyses, including pairwise 
FsT> 3D-FCA, and Bayesian clustering, suggest long-distance 
gene flow between geographically distant populations in both 
species on both coasts (Table 2; Figs. 5 and 6), all consis- 
tent with population connectivity via watercraft-mediated 
transportation of propagules. The new evidence found here, 
coupled with that from Zhan et al. (2010), suggests that in- 
traspecific genetic complexities largely resulted from genetic 
exchanges assisted by anthropogenic transfers, possibly from 
transportation of large numbers of larvae inside ballast tanks, 
and/or adults fouled on watercrafts and/or aquaculture facil- 
ities/organisms. 

Contrasting with high genetic exchange likely associated 
with human activities at both fine and large geographic scales, 
several neighboring populations exhibited significant genetic 
differentiation (Table 2; Figs. 5 and 6). Because most of the 
populations analyzed here were sampled from semienclosed 
or enclosed waters (Table 1), the high genetic differentia- 



tion between neighboring populations likely resulted from 
restricted gene flow associated with the configuration of these 
sites. Ciona larvae maybe unable to disperse naturally outside 
of semienclosed or enclosed waters. A similar genetic pattern 
has been observed in other invasive ascidians, including Styela 
clava (Dupont et al. 2009). 

Several populations exhibited lower genetic diversity and 
higher genetic differentiation in both species. For spA on the 
west coast, population MO (Monterey Bay), sampled from 
open sea shore, had lower genetic diversity than other popu- 
lations (Table 1; Table SI). Compared to other sites such as 
ports and marinas with a large number of pathways available, 
a slower natural dispersal with limited numbers of propag- 
ules can result in founder effects, leading to low genetic di- 
versity (Table 1) and high interpopulation genetic differenti- 
ation (Table 2). For spB on the east coast, all tests based on 
microsatellites including pairwise _F S t> Bayesian clustering, 
and 3D-FCA illustrated that four populations (MA, MB, 
ST, LU) were divergent from all other populations (Table 1; 
Figs. 5 and 6). Analyses of microsatellite allelic variation in- 
dicated significantly lower genetic diversity in these popu- 
lations (Table S2). The observed lower level of genetic di- 
versity and significant differentiation can result from three 
possible processes: separate introductions from genetically 
distinct sources, rapid selection due to local adaptation, 
and bottlenecks during secondary range expansion. Multi- 
ple lines of evidence suggest that the first two processes are 
unlikely to be major determinants. Detailed analysis of mi- 
crosatellite alleles did not show many unique alleles, with 
only 19 (out of 239, 7.9%) in these four populations. Re- 
gionally, given the intensive sampling and relatively low ge- 
netic differentiation of the remaining populations on the 
east coast (Table 2B; Figs. 5B and 6B), a genetically dis- 
tinct source does not likely exist on the east coast. Glob- 
ally, genetic analyses revealed a much lower level of ge- 
netic differentiation (Zhan et al. 2010). Bayesian clustering 
analysis based on a combined dataset obtained here and 
from Zhan et al. (2010) still separated these four popula- 
tions (Fig. S2). Consequently, introductions from genetically 
distinct sources both regionally and globally were unlikely. 
The lower genetic diversity at multiple microsatellite loci 
(Table SI) suggests that direct selection and/or genetic hitch- 
hiking are likely not responsible for the observed pattern. By 
contrast to the lower genetic diversity in these four popula- 
tions, one neighboring population CT (Chester marina) ex- 
hibited a level of genetic diversity similar to other populations 
on the east coast ( Fig. 2; Table 1 ) . Compared to these four sites, 
the Chester marina has both relatively higher shipping traffic 
and nearby aquaculture activities, either of which could me- 
diate gene flow. Considering the lighter commercial shipping 
to these four sites, these populations may have been acciden- 
tally seeded by limited numbers of propagules transferred 
via watercrafts, or by natural dispersal from neighboring 
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sources. Population bottlenecks during secondary range ex- 
pansions can lead to sharp decrease of genetic diversity (see 
examples in Roman and Darling 2007), and the loss of allelic 
diversity can result in a high level of genetic differentiation 
(see e.g., Zhan et al. 2009). 

Inconsistent patterns based on different 
markers 

Interestingly, we observed inconsistent patterns of genetic 
diversity when using mtDNA and nuclear microsatellites in 
both species (Table 1; Table SI). The two newly colonized 
populations of spB, BR (Brudenell River) and MR (Murray 
River), showed relatively low mtDNA diversity (n = 3, h = 
0.191), though the pattern was not repeated with microsatel- 
lite loci (Table 1 ) . In addition, the two populations, MA (Mar- 
tin's River) and ST (Stone Hurst) that possibly experienced 
a genetic bottleneck had lower microsatellite diversity but a 
similar level of mtDNA diversity to other populations (Ta- 
ble 1 ; Fig. 2 ) . A similar pattern was also detected in population 
PH (Port Hueneme) of spA (Table 1). Given the difference of 
effective population size between mtDNA and nuclear mark- 
ers, genetic drift associated with population size decrease dur- 
ing postestablishment spread may be responsible for produc- 
ing this inconsistency. In addition, different selection pres- 
sures and/or genetic hitchhiking events on mitochondrial 
versus nuclear genomes could also lead to this discordance. 

Conclusions 

Our genetic surveys add to growing evidence that complex 
genetic patterns can emerge in invasive species, including 
closely related species having similar biological characteris- 
tics. In addition, inconsistent genetic patterns also can be 
detected when using two different types of markers, that 
is, mtDNA versus nuclear microsatellites. Detailed analyses 
suggest that multiple factors, including random genetic drift 
associated with population demography change, rapid selec- 
tion due to strong local adaptation, varying propensity for 
human-mediated propagule dispersal, as well as their inter- 
actions, could be responsible for such genetic complexities. 
In addition, in order to avoid possible biased results caused 
by one type of molecular marker, our study highlights the ne- 
cessity of use of both mtDNA and nuclear markers to perform 
genetic analysis in invasive species. The complex genetic pat- 
terns observed here, coupled with multiple putative factors 
causing such complexities, open a window for us to deeply 
investigate mechanisms of rapid evolution associated with 
fast environmental changes due to frequent anthropogenic 
activities. 
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Figure SI. Bayesian inference (BI) tree (A) and haplo- 
type networks generated with TCS for Ciona intestinalis 
spA (B) and spB (C) based on a combined mitochondrial 
dataset obtained from the present study and Zhan et al. 
(2010). 

Figure S2. Bayesian inference (BI) of population structure 
for Ciona intestinalis spA (A) and spB (B) based on a com- 
bined dataset obtained from the present study and Zhan et al. 
(2010). 



Table SI. Genetic diversity at eight microsatellite loci for two 
highly invasive species, Ciona intestinalis spA (west coast) 
and spB (east coast) collected from North America. 
Table S2. P- values for the exact test of difference in allelic rich- 
ness (below diagonal) and expected heterozygosity (above di- 
agonal) using nonparameteric statistics (Mann-Whitney U 
test) for the two highly invasive species, Ciona intestinalis spA 
on the west coast (S2-1) and spB on the east coast (S2-2) of 
North America. 
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